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Abstract: In [TJ, B61 and Reese have introduced a discrete model for polymer networks 
by means of a finite element modeling. In [HIS], they present details on the method and 
a comparison with real experiments. A key parameter of their model is the size h of the 
finite element mesh, that is meant to be small in practice. The aim of the present work is to 
study the asymptotic behaviour (and the convergence of the finite element method) when 
the meshsize goes to zero. In particular, we address the properties satisfied by the model at 
the limit, depending on the properties of the mesh. 

Key-words: nonlinear elasticity, rubber, discrete system, finite element method 



* 1 DAEIMI, Universita di Cassino, via Di Biasio 43, 03043 Cassino (FR), Italy ; 2 Dipartimento di 
Matematica e Applicazioni 'R. Caccioppoli', Universita di Napoli, via Cintia, 80126 Napoli, Italy ; 3 CER- 
MICS - ENPC & INRIA Paris- Rocquencourt, 6 et 8 avenue Blaise Pascal - Champs sur Marne, France ; 
Contact: alicandr@unicas.it, cicalese@unina.it, antoine.gloria@inria.fr 



R. Alicandro 1 , M. Cicalese 2 , A. Gloria 3 B 



Theme NUM — Systemes numeriques 
Projet MICMAC 



Rapport de recherche n° 6270 — August, 2007 — [9] pages 



Unite de recherche INRIA Rocquencourt 



Proprietes de convergence du modele discret de 
Bol-Reese pour le caoutchouc 

Resume : Dans p], B61 et Reese ont introduit un modele discret de reseau de polymeres 
par elements finis. Dans [21 Sj, les auteurs detaillent leur methode et comparent les resultats 
numeriques obtenus a des experiences reelles. Un des parametres cles de leur modele est la 
taille du pas h du maillage, qui est petit en pratique. L'objectif de ce travail est d'etudier 
le comportement asymptotique (ainsi que la convergence de la methode des elements finis) 
lorsque h tend vers zero. En particulier, nous etudions les proprietes satisfaites par le modele 
a la limite, en fonction, des proprietes du maillage. 

Mots-cles : elasticity non lineaire, caoutchouc, systeme discret, methode des elements 
finis 
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Phenomenological constitutive laws for rubber-like materials often involve parameters 
which lack of physical motivation. In addition, they are usually difficult to fit in practice. 
Much attention has been paid in the recent years to microscopically-based models. This 
kind of models aims at pointing out the microscopic features that govern the macroscopic 
behavior of the material, such as the geometry of the underlying polymer network. Such a 
model has been introduced by B61 and Reese in pp. In addition to the physical parameters 
provided by physics at the scale of the polymer chains (such as the free energy of a chain) , 
this model involves two geometric parameters: a typical distance h > (which is meant to 
be small) and a tetraedral mesh. 

The aim of the present work is to study the asymptotic properties of the model when h 
goes to zero for different assumptions on the mesh. 

This note is organized as follows: 
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1 The B6L-Reese model 
1.1 The geometry 

For the sake of simplicity, let denote the unit cube of M 3 and Tj, be a tetraedral mesh of 
O. The typical distance we associate with the mesh is defined as follows: 




where iV j denotes the number of elements (tetraedras) of % L . 

From a modeling point of view, each edge of the mesh represents an elastic spring with 
an energy related to the underlying polymer network. 

Let d be a displacement field of the mesh defined at each vertex Xi of 7^ by c£j. One can 
associate with d a continuous displacement field on SI as follows: Vh is the unique piecewise 
affine function on 7^ such that for each vertex ir, of 7^, Vh(xi) = Xi + rf,. 

The energy of the system deformed by Vh is then given by 

Ten (*„*,)eT v lXl J Ja 

where T denotes elements of Th, (xi,Xj) are distinct pairs of vertices in the tetraedra T, 
Wfj is a pair potential, is a scaling factor and W vo i is an energy density accounting for 
volume changes of the deformed system. 

Let us now describe the pair-potential energy. 



1.2 The energy 

The energy of the elastic springs is related to the energy of a polymer chain as follows: 



+ ln- 



sinh£ 



c 

0' 



(2) 



where j3 is the inverse of the absolute temperature, £ 1 is the inverse of the Langevin 
function, k and c are constants and nfj is a typical number of segments of a polymer chain, 
can be interpreted as the free energy of a polymer chain made of nfj segments at the 

length \J n Jjl, where I denotes the length of a segment. 

The factor ff 3 is a measure of the number of polymer chains per unit of volume. We will 
typically consider it as constant on fl. 



INRIA 



Convergence of a discrete model for rubber 



5 



Finally, the energy accounting for volume changes is given by 

W vo i(0 = * (J 2 -l-lnJ), (3) 

where K > 0, £ is a deformation gradient and J = det(£). The physical origin of this 
contribution is the Van der Waals forces, which prevent atoms from being too close to one 
another. The scale of this forces is far larger than the scale of the polymer network, which 
explains why this contribution is already coarse-grained and somewhat uncorrelated with 
the description of the network. 

In the following section, we study the convergence when h — > of the sequence of 
minimization problems 

mi{E h {v h ),v h G V h +BC}, 

where V% denotes the space of piecewise afHne functions on Th and BC stands for the 
boundary conditions (let say mixed imposed displacement and free traction for instance). 
In particular, we address both the convergence of the infimum of the energy and the con- 
vergence of the minimizers (the deformation field) when h — ► 0. This convergence analysis 
is performed using T-convergence (see [11 [5] e.g.). 



2 Convergence analysis 

2.1 Assumptions on the energies 

For technical reasons, in what follows, we will assume that: 

• Wy satisfies the following standard growth condition of order p > 1 for all T and ij: 
there exist C > c > such that for all r, 

c\r\* - 1 < W%(r) < C(\r\> + 1), (4) 



• W vo i satisfies Q from above. 

In particular, one cannot directly deal with the inverse of the Langevin function. To 
satisfy the growth condition ([4]). one may use a truncated series expansion of C~ l , as it is 
done in [DEI [3]: 

n 1/ N n 9 3 297 5 1539 7 9n 

c-\ p) = z P+T y + -? + — p T + o {p »). 

Replacing £ _1 by its development up to order 7 in (El, we obtain a free energy density for 
a polymer chain which satisfies (J4j) for p = 8. 

The energy density W vo i given by (J3j) does not satisfy |4]) since lim dGt (£)^ W vo i(£,) = +oo 
while det(£) — > does not imply |£| — * oo. One possible solution consists in taking a cut-off, 
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which amounts to relaxing the constraint of incompressibility. A typical cut-off reads as 
follows: 

W*> (ft - / Ww *@ if det( ^ >V K\ 

WvolU) - | K( v 2_ 1 _ lnr?) otherwise. [ °> 

For 77 > fixed, W^ ol satisfies (j4]) for C v large enough. To recover the incompressibility 
behavior, one can first perform the convergence analysis on the discrete to continuum process 
for 77 > and then let r\ go to zero at the continuum level. We refer the reader to [51 [7] for 
technical details. 

We are now in position to state our convergence results. 
2.2 The periodic case 

In this paragraph, we address the convergence of the Bol-Reese model in the periodic case. 
To this end, let us assume that the edges of Th are obtained by the periodic replication of 
a unit cell of edges (up to border effects on the boundary) , as illustrated in two dimensions 
on Figure [TJ 



Figure 1: Two-dimensional periodic mesh, the unit cell is in bold 

Let us further assume that the energy of the elastic springs does not depend of the unit 
cell and on the mesh lengthscale h. As a particular case, one can take nfj to be constant 
on Th, as it is done in [HEIIH]. Then, [HI Theorem 4.1] shows there exists an homogeneous 
quasiconvex energy density Whom, which is frame- invariant and satisfies a growth condition 
of order p, such that: 

• lim inf {E h (v h ),v h <E V h + BC} = inf {E hom (v),v G V + BC}, where V is the Sobolev 
space W^ P (Q), and 

E hom (v) = / W hom (Wv); 
Jn 

• Whom is given by the asymptotic discrete homogenization formula 

W hom {C) = lim mf{E h (v h ),v h G V h ,v(xi) =£-Xi ifd(x u dn) < 2h}. 
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In addition, minimizers u% of Eh on Vh + BC weakly converge in W 1:P (il) to minimizers 
Uhom of Ehom on V + BC . This proves in particular the convergence of the solutions of the 
finite element modelling of rubber introduced by B61 and Reese. 

In the periodic case, the energy density at the limit is frame-invariant. However it is not 
isotropic in general. To ensure the isotropy, one may use a stochastic framework. The end 
of this paragraph is dedicated to the study of a simple example for which the energy density 
at the limit is not isotropic. 

Let us consider simple linear springs by setting 



on the mesh sketched on Figure[TJ Then, recalling [HI Remark 5.2], the homogenized energy 
density is quadratic and we end up with a convex minimization problem on one single 
periodic cell. The infimum is trivial and shows that the homogenized energy is not isotropic. 
It is actually enough to see that the material is stiffer in the direction e-2. — ei than in the 
direction e\ + e-i (where e\ and e2 denote the canonical basis of M 2 ). 

2.3 The stochastic case 

Let us quickly recall the concept of admissible stochastic network on a probability space 
(S, !F, P) used in [HE], which is a particular case of the stochastic networks introduced by 
Blanc, Le Bris and Lions in [TQlfTT] . 

Let A = {yi}i e z d £ (K^) 2 be a set of points. We say that A is an admissible set of 
points if it satisfies the two following conditions: 

i. there exists R > such that #A n B(y, R) > for all y e R d ; 

ii. there exists r > such that d(yi, A \ {yi}) > r for all i £ Z d . 

In particular, to each admissible set of points A one can associate a Delaunay triangulation 

T>(A), A stochastic lattice L : 3 — > (M d ) Z is said to be admissible, if for P-almost every 
lo € 5, C(ui) is an admissible set of points, and if the Delaunay triangulation is regular in 
the sense of the interpolation theory. 

Given a realization A(w) of a stochastic lattice, one may rescale the lattice by a factor h, 
setting = hyi(u). Let us then define A(uj) h (tt) = {y' t l e R 3 ,yi € A(u>) and y£ € fi}, 

which is the intersection of the rescaled lattice with ft. Using the notations of B61 and Reese, 
one may think of a mesh Th as the Delaunay triangulation of A(oj) h (il). 

In what follows, we assume that Th is related to a stochastic network in the way described 
above and we add an index u> to make the stochastic dependence more explicit. Futhermore 
we will make some hypotheses related to the stochastic network and the probability space: 
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(a) there exists a stationary translation group that acts on the stochastic lattice and which 
is ergodic for P; 

(b) there exists a stationary group of rotations that acts on the stochastic lattice (the 
stochastic lattice is then said to be rotation invariant on average). 

Roughly speaking, the translation invariance allows us to obtain a deterministic limit whereas 
the invariance by rotation implies that the energy density at the limit is isotropic. We refer 
the interested reader to [7] for the precise formulation of the assumptions and the proofs of 
the following result. 

Within hypothesis (a), there exists an homogeneous quasiconvex energy density Whom, 
which is frame-invariant and satisfies a growth condition of order p, such that: 

• For almost all lo e S, 



• Whom is given by the asymptotic discrete homogenization formula 

Whom(0 = lim/wo Js inf < £ V£,v(x t ) = £ ■ x, 

tfd(xi,d£l) < 2hR}dP{uj). 

In addition, minimizers of E% on Vff + BC weakly converge in W 1,p (£l) to minimizers 
Uhom of Ehom on V + BC. This proves in particular the convergence of the solutions of the 
finite element modelling of rubber introduced by B61 and Reese. 

In addition, if (b) holds, the stochastic network is rotation invariant on average, and 
Whom is isotropic. 

3 Conclusion 

In the first section, we have recalled the finite element modelling of rubber developed by B61 
and Reese and we have pointed out some convergence issues related to their formulation. 
In the following sections, we have analyzed the asymptotic behaviour of the model. In 
particular, we have identified two cases for which their model provides a finite element 
approximation of an underlying continuous model. The continuous models obtained both 
satisfy the frame invariance property and the usual incompressibility behaviour used for 
rubber modelling (up to some technicalities detailed in [2]). However, only the stochastic 
model is ensured to yield an isotropic energy density at the limit. We have also provided a 
periodic example which is not isotropic at the limit. 



lim inf v% G V? + BC} = inf {E hom (v), veV + BC}, 

h— >0 

where V is the Sobolev space W 1,p (il), and 
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